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Abstract 

The production and space-time evolution of charm and bottom quarks in nucleus-nucleus colli- 
sions at RHIC and LHC are investigated with the partonic transport model BAMPS (Boltzmann 
Approach of MultiParton Scatterings). Heavy quarks, produced in primary hard parton scatterings 
during nucleon-nucleon collisions, are sampled using the Monte Carlo event generator PYTHIA 
or the leading order mini-jet model in conjunction with the Glauber model, revealing a strong 
sensitivity on the parton distribution functions, scales, and heavy quark mass. In a comprehensive 
study exploring different charm masses, K factors, and possible initial gluon conditions, secondary 
production and the evolution of heavy quarks are examined within a fully dynamic BAMPS sim- 
ulation for central heavy ion collisions at RHIC and LHC. Although charm production in the 
quark- gluon plasma can be neglected at RHIC, it is significant at LHC but very sensitive to the 
initial conditions and the charm mass. Bottom production in the quark-gluon plasma, however, is 
negligible both at RHIC and LHC. 
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I. INTRODUCTION 



Experimental observations from relativistic heavy ion collisions at the BNL Relativistic 
Heavy Ion Collider (RHIC) indicate the production of a hot and dense partonic medium 
MiM, commonly referred to as the guark-gluon plasma (QGP). A good agreement between 
ideal hydrodynamic simulations [5 T ] and measurements js), [J of the elliptic flow hints at 
fast thermalization of this system, which behaves like a nearly perfect fluid, that is, possesses 
a very small shear viscosity to entropy density ratio jiol . 111 ]. 

Heavy quarks, that is, charm and bottom, are a unique probe for this medium. Due to 
their large mass (Mq 3> A QCD , Q = c,b), a large amount of energy is needed to produce 
heavy quarks. Such high energy densities are primarily found at the early stage of heavy 
ion collisions: in hard scatterings of partons in the nucleons of the heavy ions or during 
the early phase of the QGP. In that ener gy d omain, the running coupling of the strong 



12j. Therefore, heavy quark production should 
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14j . even for small 



interaction a s is small and nearly constant 

be describable within the framework of perturbative QCD (pQCD) 
transverse momenta. 

Another theoretically proposed implication of the large mass is the "dead cone effect" 
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16], which implies a smaller energy loss via gluon radiation compared to light quarks 



and delays the thermalization of heavy quarks by a factor of ~ Mq/T 



171 ] . resulting in 



a thermalization time o: 



observed energy loss jl8, Q and elliptic flow 



about the lifetime of the QGP. In contrast, the experimentally 



19] of heavy flavor is comparable to that 
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25] . Due to the 



of light quarks, the reason for this puzzle being under investigation 
- in principle - unique identification of heavy quarks because of their flavor and roughly 
known production time during the early stage of the collision, their distributions can reveal 
information about the interaction history, rendering them as an ideal probe of the medium. 

In principle, heavy quarks can be produced at three stages of the collision: during initial 
hard parton scatterings, in the QGP, or during the hadronic phase. As we will show in 
this article, most of the heavy quarks are created in hard parton scatterings during initial 
nucleon-nucleon collisions in the heavy ion collision at RHIC. However, at the CERN Large 
Hadron Collider (LHC) secondary production during the QGP phase becomes important. 
Here, one often distinguishes between prethermal and thermal production depending on 
whether the partonic medium is already thermalized or not. In the present article, we focus 



on the first two possibilities of heavy quark production and neglect the hadronic phase, 
which hardly contributes to the heavy quark yield. 

An interesting signature in heavy ion collisions is the suppression 26] or enhancement 



of heavy quarkonia like J/tp (hidden charm) [27]]. In this article, however, we examine only 
open heavy quark production and postpone the investigation of hidden charm or bottom 
to future investigations. For the simulation of the heavy quark production in the QGP, we 
use the partonic transport model called the Boltzmann Approach of MultiParton Scatterings 
BAMPS) and study the impact of various initial gluon distributions such as from PYTHIA 



28j |. the mini-jet model 29|, |30(, and the color glass condensate 3l|, |32j. The initial heavy 



flavor yield in hard parton scattering is obtained with PYTHIA and compared with that 
from leading order (LO) pQCD. 

Most of the time we will only refer to charm quarks instead of heavy quarks in general. 
However, the majority of the concepts and findings apply qualitatively also to bottom quarks 
if one takes the larger mass into account. For LHC calculations, we will explicitly mention 
bottom production and make predictions on the bottom yield. 

This article is organized as follows. First we introduce our model BAMPS and the par- 
tonic cross sections for heavy quark production. In Section UTTl we use this model to estimate 
the chemical equilibration time of a static system filled with gluons and charm quarks and 
compare it to the analytic solution. The next two sections address charm and bottom produc- 
tion in heavy ion collisions. After investigating the initial heavy quark distributions from 
primary nucleon-nucleon scatterings in Section \W\ we study the heavy quark production 
during the QGP phase in Section [V] Here, the charm quark fugacity is also calculated and 
related to the estimated chemical equilibration time from Section II III Finally, we conclude 
with a short summary in Section IVII 



II. PARTON CASCADE BAMPS 



For the simulation of the QGP phase we use the partonic transport model BAMPS 33, 



34|, which stands for Boltzmann Approach of MultiParton Scatterings. BAMPS simulates 



the fully 3 + 1 space-time evolution of the QGP at heavy ion collisions by solving the 



Boltzmann equation, 



— + PlA ] f.( r D t ) = C 2 ^ 2 + C 2 ^ 3 4 



(1) 



dynamically for on-shell partons with a stochastic transport algorithm and pQCD inter- 
actions. Cj are the relevant collision integrals, and /j(r, p i; t) the one-particle distribution 
function of species i = g, c, c, since light quarks are not included yet. Bottom quarks can 
also be investigated if one changes c — > b and M c — > M^. In addition to the binary collisions 
2—7-2, also 2 -H- 3 scatterings for the gluons are possible. That is, the following processes 
are implemented in BAMPS: 



9+9^9+9 
9+9^9+9+9 
9+9+9^9+9 
9 + 9 -> c + c 
c + c-> g + g 
9 + c-t g + c 
9 + c^f g + c 



(2) 



Within the stochastic method, the probability for a collision of two particles during a time 
step At in a volume element AV can be obtained from the collision term of the Boltzmann 
equation [331 ] : 

At 



AV 



(3) 



where 



Vre\ 



^(PrP 2 ,) 2 - m( 
E\E 2 



2 mi 



(4) 



stands for the relative velocity, and <y 2 ^2 for the binary cross sections. The cross section for 
99 ~~ 99 as wen as the treatment of the gg O ggg interactions are given in [33 ] . 



The most dominant charm production process is gluon fusion g + g — > c + c with a 
differential cross section of 



dcr, 



gg-^cc 



\M 



gg^cc\ 



dt 



167TS 2 



(5) 



\Ai gg ^,cc\ is the averaged matrix element, that is, the averaged (summed) over color and 
spin of the incomin g (o utgoing) particles. It can be expressed in terms of the Mandelstam 



variables s, t and u 



36 



37j: 



\M 



gg^cc\ 
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7r 2 a 2 



(M 2 -t)(M 2 -u) + 



8 (M 2 -u M 2 — t 
+ 



16M 2 



M 2 + t 



+ 



3 \M 2 -t 
M 2 + u 



M 2 



1 1 



_(M 2 -t) 2 (M 2 -u) 



-(2M 2 

s 



u 



+ 



6 M 2 (t-u) 



2 M 2 (s-AM 2 



s (M 2 - t)(M 2 - u) 3(M 2 -t)(M 2 -u) ' ^ 
M denotes the mass of the charm (or bottom) quarks. If not otherwise specified, we use 



Mr 



1.5 GeV for charm and Mh = 4.8 GeV for bottom quarks, which are widely adopted 



in the literature 



13, 



27 



36 



38 



44j . For the coupling of the strong interaction we choose to 
take a constant value of a s = 0.3. 

After substituting u with the relation s + t + u = 2M 2 , the total cross section is obtained 



by integrating Equation (5) , [37J, |45M48| 



3s 



AM 2 M 4 
1 + + 



log 



1 + X 



31M 5 
4s 



X 



with the abbreviation 



X 



4M 2 



(7) 



(8) 



The angular distribution of the charm and anti-charm quarks after the collision is sampled 



using Equation (5) 



The back reaction c + c— > g + g is negligible in heavy ion collisions at RHIC and LHC 



271 ] due to the small number of produced charm quarks compared to the number of gluons. 



However, we want to implement the back reaction in a box model in Section II III in order 
to estimate the chemical equilibration time scale. The cross section of this process can be 
obtained through detailed balance, 

1 64 1 



®cc^gg ~ 2 9 -)/ 2 °"^ CC 



(9) 



The factor 1/2 comes into play due to gluons being identical particles. 64/9 takes the color 
and spin averaging and summation into account whereas 1/x 2 is a kinematical factor. 

The process for charm production by light quark and anti-quark annihilation q + q — > c+c 
is not yet included in the cascade. However, we use this process to estimate the initial charm 



quark yield in heavy ion collisions (see Section TlV Cp . The cross section for that process is 
given, for instance, in 37], |45|-|48|. 



III. CHARM PRODUCTION IN A STATIC MEDIUM 

In this section, we test our implementation with BAMPS by comparing our numerical 
results to the analytic solution of a rate equation. To keep the problem simple and find an 
analytic solution, we consider a box of gluons and charm quarks, in which only charm-anti- 
charm production through gluon fusion and the back reaction 



9 + 9 -> c + c 
c + c-> g + g 



(10) 



are allowed, whereas all other interactions between gluons and charm quarks are forbidden. 
However, we checked that adding other possible processes like gg —> gg, ggg -H- gg, gc —> gc, 
and gc — > gc does not have an impact on our findings. The additional 2—7-2 processes 
ensure the kinetic equilibration and ggg -H- gg maintains the chemical equilibration of the 
gluons, but they do not affect the chemical equilibration time scale of the charm quarks. 



evolution of the charm density 



36 



48 



450| 



For the system with only processes from (10), we write down a rate equation for the 



<9 M (n c w M ) = R, 



gg-^cc 



R 



cc-^gg i 



(11) 



where R gg ^ C c and Rcc^gg denote the rates and = 7 (1, v) (7 being the Lorentz factor) 
the four-velocity of the considered volume element with charm density n c , which is identical 
to the charm pair density n c5 since charm and anti-charm quarks are always produced in 
pairs. 



The rates are given by 



36 



49 



50| 



R 



gg^cc 



- o- 



gg^cc 



v T ei) n. 



2 1 2 

" : 2^ 



R, 



cc^gg 



cc-^gg 



Vre\) n c n- c 



(12) 
(13) 



In Equation (12), n g stands for the gluon density and the factor 1/2 is needed to take into 
account that gluons are identical particles. a g := (a gg ^ C c v re {) and a c := {a C c-> gg v TC \) are the 



mean cross sections weighted with the relative velocity t> re i defined in Equation (4) In all 



calculations we use a constant coupling of a s = 0.3 for the cross sections, which are given 
in equations (7) and (9) 



A. Analytic solution of the rate equation 



For a static box the four-velocity of each volume element is = (1, 0, 0, 0) and the rate 



equation (11) simplifies to 



dt'lT'cc Rgg-^tcc Rcc-^gg <^9^g ^c^cc 



(14) 



In chemical equilibrium the rates are equal, R e g ^ cS = Rl^ gg . From that relation one can 



obtain the charm density in chemical equilibrium, 

1 "-tot 



n. 



eq 



+ 1 



(15) 



yj2<Tg/(T C 

where n tot = n g + 2n c5 denotes the constant total particle density. 

Initially the box may only contain gluons which are chemically and thermally equilibrated. 
That is, the initial gluon density for an initial temperature of T reads 



%(* = 0) = ^ 



(16) 



(gluons are treated as Boltzmann particles, v g = 2 ■ 8 = 16 is the degeneracy factor for 
gluons), whereas for charm quarks n c5 (t = 0) = 0. 



Solving the rate equation (11) provides the time evolution of the charm density, from 



which one can read off the chemical equilibration time scale. Taking into account that the 
total particle number is constant, the charm quark density as a function of time is given by 



/ \ 1 ^tot 
n cc{t) = ~- 



2 1 -C 



Jt/r 



(C + i) -c + i 



e 2 */- U + 1 



(17) 



where the abbreviations 



n t q "tot - 2n 



eq 

cc 



2n. 



eq 



2n, 



cq 



2ri 



cq 



2n 



cq 



crgntotn e g q 



2n tot n, 



eq\ 



(18) 
(19) 



have been introduced. 
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The solution is implicitly dependent on the temperature via a c and a g and is only valid 
assuming that the temperature stays constant over the whole period of time. As we will 
see in the next section, due to the mass creation this is, however, not exactly but in a good 
approximation the case. 



B. Comparison between analytic and numerical solution 



The numerical solution is obtained with the parton cascade BAMPS (see Section [IT])- In 
order to compare it to the analytic solution from the previous section only the two processes 



from (10) are allowed. Again, the initial charm density is n c5 (t = 0) = and the gluons are 



sampled thermally employing n g (t = 0) according to Equation (16) with temperatures of 
400 MeV and 800 Me V, which correspond to the expected temperatures of the quark-gluon 



plasma at RHIC |5lM57j and LHC [53|, |56|, |57| , respectively (cf. also Figures fT9l and l2Tj) . 



1. Charm production at RHIC temperature 



The initial temperature of the gluons is taken as T g = 400 MeV [5lH57j . Being precisely, 
during the evolution of the box the temperature is ill-defined because of the gluons being not 
in thermal equilibrium (due to the neglect of elastic scatterings among the gluons). However, 
one can extend the temperature definition to this non-equilibrium domain, introducing an 
"effective gluon temperature" 

EJt) e (t) . , 

Here E g (N g ) stands for the total gluon energy (number) in the box, and e g and n g for the 
respective densities. 



The time evolution of the gluon temperature is illustrated in Figure 1 The temperature 
declines to about T g = 370 MeV since part of the kinetic energy is used for the production 
of massive charm quarks. Consequently, this implies a problem in comparing the numerical 
results with the analytic solution, which is only valid for a constant temperature. Therefore, 
we will compare the analytic solutions for initial and final temperatures against the numerical 
results at the beginning and at the end, respectively. 



Figure 2 shows the time evolution of the charm pair density n c£ . In addition, the analytic 
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Figure 1. (Color online) Time evolution of the gluon temperature in a static system with an initial 
temperature of To = 400 MeV. The temperature decreases due to the creation of massive charm 
quarks. 



solutions of the rate equation for the charm density from Equation (17) for initial and final 
temperature are also plotted. In the beginning, the numerical results are in very good 
agreement with the analytic solution for T = 400 MeV. Thereafter, the temperature of the 



gluon plasma decreases (cf. Figure 1 ), and the numerical results must be compared with the 
analytic solution for the stationary final temperature of T = 370 MeV, both being in good 
agreement. 



The rates for gg — > cc and the back reaction are shown in |Figure 3} As one would expect 
from the principle of detailed balance, the rates converge to a mutual constant equilibrium 
value with time. The chemical equilibration time scale may be defined as r cq . By fitting the 
curve m 



Figure 2 with 



n C c{t) 



-t/T e , 



(21) 



we obtained an approximate value of 700 fm/c for r eq . Of course, this result of the chemical 
equilibration time was estimated within a simple static model without expansion. However, 
the order of this result should be roughly the same as in realistic heavy ion collisions, which 
indicates that not very many charm quarks are produced during the QGP phase at RHIC 
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Figure 2. (Color online) Time evolution of the charm pair density n c g in a static medium with an 
initial temperature of To = 400 MeV. For comparison, the analytic solutions from |(17)| for initial 
and final temperature are also shown. 
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Figure 3. (Color online) Time evolution of the rates in a static system with an initial temperature 
of Tq = 400 MeV. Again, the analytic solutions for initial and final temperatures are also plotted. 
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Figure 4. Temperature dependence of the chemical equilibration time T eq) T =cons t. for charm pro- 
duction in a static medium of a constant temperature. 

considering this large time scale of chemical equilibration compared to the lifetime of the 
QGP, as we will explicitly show in Section IV Al 

An obvious expression for the temperature dependence of be obtained by differ- 



entiating Equation (21) by t, equating with Equation (14), and evaluating at t = 0: 



cq 



R, 



gg^rcc\t=0 



(22) 



It is related to r from Equation (19), which does not have a direct physical meaning, by 

1 



' eq 



T 



(23) 



1 + ^2a g /a c 

The temperature dependence of r eq can be obtained analytically with the assumption of a 
constant temperature throughout the evolution of the system. In order to distinguish it 
from the r eq obtained numerically with BAMPS, in which the temperature drop due to the 
mass creation is considered, we label the analytic solution as r eqi T= C onst. 



Figure 4 depicts the 



temperature dependence of r eqi T=const. • Compared to the results from BAMPS, these values 
for a constant temperature are slightly larger. For high temperatures, which can be present 
at LHC at a very early stage of the QGP phase (cf. |Figure 21 ), the chemical equilibration 
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Figure 5. (Color online) Evolution of gluon and charm quark fugacities in a static system. The 
final value is above 1 due to the fixed number of particles and the temperature drop. 

time lies in the same range as the lifetime of the QGP. Consequently, a substantial charm 
production in the QGP at LHC can be expected, if the initial temperature of the medium 
is large. We will address this in more detail in Section IV Bl 

The fugacity, which is defined by 



A<(t) 



nj{t) 

chem.eq 



(24) 



for particle species i, is an interesting variable for investigating the chemical equilibration. 



Figure 5 shows the fugacities of gluons and charm quarks in our simulation as a function 



of time. Due to detailed balance, they adopt the same value in equilibrium, although that 
differs from 1, a consequence of the total particle number being constant. As a note, explicitly 
allowing inelastic processes (gg -H- ggg) as in the full transport simulation (see Section |V|) 
ensures chemical equilibration and leads to a final fugacity of 1 both for gluons and charm 
quarks. 
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Figure 6. (Color online) Evolution of the gluon temperature in a static medium with an initial 
temperature of T = 800 MeV. 



2. Charm production at LHC temperature 



For the initial gluon temperature, a value of T„ = 800 MeV 



53|,|56|,|57j (cf. also|Figure 2l) 



is chosen, which drops to about T„ = 720 MeV as shown in Figure 6 due to the mass creation 



Figure 7 depicts the numerical results of the charm pair density evolution, which is in 



excellent agreement with the analytic solutions taking the temperature drop into account. 



From that figure, using Equation (21), the time scale of chemical equilibration r eq can be 
estimated to about 30 fm/c, which is much smaller than for RHIC temperature, but still 
sizable. 



The rates for both considered processes can be found in Figure 8 



IV. INITIAL PARTON DISTRIBUTION IN HEAVY ION COLLISIONS 

The initial distributions of the partons play a crucial role for the dynamics of the heavy ion 
collision. In this section, we want to outline a prescription to describe the heavy ion collision 
- according to the Glauber model - as a superposition of nucleon-nucleon collisions, which 
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Figure 7. (Color online) Charm quark density as a function of time in a static medium with an 
initial temperature of To = 800 MeV. The analytic solutions for initial and final temperature are 
also shown. 
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Figure 8. (Color online) Evolution of the rates in a static medium with Tq = 800 MeV. 
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are sampled with the event generator PYTHIA 28] . In that framework, we also study charm 
production and compare the yields with experimental data from nucleon-nucleon collisions. 
To be able to examine the impact of the initial conditions on charm production during the 
QGP phase in Section |V] we also discuss other models for the initial parton distributions 
such as the mini-jet model and the color glass condensate. 

The prescription for the position sampling of the partons according to a geometric model 



is described in great detail in 



33|. 



A. Parton and momenta sampling with PYTHIA 



For our simulation we use PYTHIA 6.4 28|], allow hard and soft QCD interactions, and 



turn off the hadronization. PYTHIA distinguishes between soft and hard events. Therefore, 
we define particles stemming from a hard (soft) event as being hard (soft), regardless of 
their momenta. The only exception to this rule is that beam remnants such as diquarks are 
always considered as soft. 

For nucleon-nucleon collisions at RHIC with a center of mass energy of ^/stvtv = 200 GeV, 
PYTHIA yields the following results. On average, 53% of all processes are hard and 47% 
soft. All particles produced in hard processes are partons. These hard partons account 
for about 58 % of all particles created in nucleon-nucleon collisions. On the contrary, in 
soft processes, most formed particles are non-partonic, for instance, diquarks or excited 
nucleons. Their fraction is about 28 % of all particles, whereas partons which are produced 
in soft events only account for 14 % of all particles. The total number of produced particles 
in one nucleon-nucleon collision averages to 7.5, whereof 5.4 are partons. The total energy 
is, of course, ^/snn = 200 GeV, but all partons together possess only 30 GeV. Consequently, 
non-partonic particles account for 85 % of the total energy. The reason for that is their large 
energy per particle ratio, which results from the non-partonic particles either being beam 
remnants or stemming from elastic or diffractive hadronic processes and, therefore, carrying 
a huge amount of energy. 

The reason for the separation in soft and hard particles lies in the different scaling be- 
havior from nucleon-nucleon to heavy ion collisions. According to the Glauber model, hard 



processes scale with the number of binary collisions, 58 1 



A bin (b) =a p+p T AB (b) , (25) 
15 



where b is the impact parameter, T^^b) the nuclear overlap function for collisions of two 



nuclei A and B, and cr p+p ~ 40 mb for RHIC and cx p+p ~ 60 mb for LHC the ine 



astic p+p 



cross sections. For central collisions, we approximate Taa{^> = 0) = A 2 /tcR\ 3jJ |59j, which 
leads to Taa(^ = 0) = 30.4 mb -1 for Au+Au collisions at RHIC being in excellent agreement 
with a numerical calculation using [60]. Therefore, 

iV bin = a p+p T AA (b = 0) « 1200 (26) 

is the number of binary collisions at RHIC for ^/sJjn = 200 GeV and b = 0. Taking 



shadowing into account reduces this number to about iVbin ~ 1000 [6JJ, which we use as 
a scaling factor for hard partons, Ch ar d := N^ in . To obtain the scaling factor C so f t for soft 
particles, one can make use of energy conservation, 

EcM = -Epp.hard Chard + -Epp, soft C so ft , (27) 

where Ecu stands for the total energy available in the heavy ion collision (at RHIC Ecu = 
200 A GeV = 39 400 GeV) and E pP) hard(soft) for the total energy of all hard (soft) particles in 
one nucleon-nucleon collision. Solving this equation for C so f t leads to C so ft ~ 100, which is in 
the same order as the number of nucleons, the usually used scaling factor for soft particles. 



Employing the scaling prescription obtained from Equation (27) results in the following 
particle yields for Au+Au collisions at RHIC with ^s^n — 200 GeV: Partons from hard 
processes account on average for about 93 % of the produced particles. 4.5 % of all particles 
are non-partonic particles from soft processes, whereas 2.5 % are soft partons. The total 
number of produced particles is about 4600, of which 4400 are partons. However, for the 
energy distribution the picture looks a bit different. Hard partons carry just 53 % of the 
total energy and non-partonic particles with 44 % nearly the same amount. Soft partons only 
account for 3 % of the total energy. Consequently, the whole energy deposited by partons, 
which is available for the parton cascade, is about 56 % of the total energy. 



Figure 9 depicts the rapidity distribution of the particle number and of their transverse 
energy in a Au+Au collisions at RHIC. Gluons from hard processes dominate the spectrum 
at mid-rapidity in both distributions. In addition, hard quarks take a considerable fraction 
of the transverse energy. Partons from soft processes have lost much of their significance 
compared to unsealed nucleon-nucleon collisions due to the smaller scaling factor of soft 
particles. 
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Figure 9. (Color online) Rapidity distribution of particle number dN/dy (a) and transverse energy 
dEx/dy (b) in central Au+Au collisions at RHIC. As explained in the text, the distributions are 
obtained by simulating nucleon-nucleon collisions with PYTHIA and scaling them to heavy ion 
collisions. CTEQ61 is used for the parton distribution functions, "others" denotes non-partonic 
particles such as diquarks or protons from soft processes. The distribution of soft gluons is almost 
zero and is, therefore, not shown here or in the following figures. 
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Figure 10. (Color online) Transverse momentum spectra dN / {pxdpxdy) of particles at rapidity 
y G [—0.5,0.5] in central Au+Au collisions at RHIC. 



In Figure 10 the transverse momentum spectra are shown. Of course, hard partons also 
dominate these spectra. The transverse momentum of soft quarks is always smaller than 
2 GeV because of a momentum cut-off in PYTHIA. For partons from hard events, we did 
not introduce a cut-off. However, PYTHIA has an internal cut-off for so called semi-hard 
scatterings, which we also consider as hard events. The cut-off is dependent on y/s and lies 
for RHIC energy at 141 MeV causing a jump at this value in the pt spectra of hard partons. 

For LHC, the initial parton distribution is also sampled with PYTHIA and scaled by 
using the same prescription. For central Pb+Pb collisions at y / s/v r /v = 5.5 TeV, the overlap 
function is TaaO* = 0) ~ 32.7 mb -1 and proton-proton cross section c p+p rj 60 mb, which 
determines the number of binary collisions to be about 2000. Actually, we reduce this value 
to A^in ~ 1500 to take also shadowing into account |62|, |63|. |Figure lT] depicts the rapidity 
distributions of the particle number and their transverse energy for heavy ion collisions at 
LHC. 

Currently, BAMPS does not include light quarks. Therefore, to conserve energy and 
particle number, light quarks from PYTHIA are converted to gluons before we evolve the 
system with the parton cascade. This conversion does not account for effects such as the 
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Figure 11. (Color online) As in Figure 9 but for central Pb+Pb collisions at LHC (V-sjvn 
5.5 TeV). 



faster thermalization of gluons compared to light quarks and the larger cross section of 
gg — > cc than of qq — > cc. However, since gluons dominate the early phase of the QGP, 
where most of the heavy quark production takes place, we do not expect an impact on the 
production of heavy quarks. Nevertheless, this will be checked in the near future after light 
quarks are included in BAMPS. 
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B. Other models for the initial parton distribution 



Two other models for the initial par ton distribution were already used in the parton 



model 



64 



65) 



cascade BAMPS: the mini-jet model [33|, |34j and a color glass condensate (CGC) inspired 



In the mini-jet model 29|, |30), the initial parton distribution is given by several indepen- 



dent 2-jet events, which are sampled according to |66j 



da 



jet 



dp%dyidy 2 



K^2x 1 f a (x li PT)x 2 fb(x2,PT) 



da. 



ah 



a.b 



dt 



(28) 



where pt denotes the transverse momentum, y the rapidity and x the Bjorken variable. The 
cross section <r ab is calculated in LO pQCD and a K = 2 factor is introduced to account for 
higher orders. 

To avoid problems at low momenta where pQCD is not valid anymore, we employ a 
momentum cut-off, which is set to po = 1.4 GeV for RHIC in order to get the final transverse 



energy distribution of the gluons in agreement with data 
3.3 GeV to get a gluon yield, which is comparable to the PYTHIA result. 
The number of produced partons is given by 

N part0 ns(b) = a ict T AB (b) . 



67]]. For LHC we choose po 



(29) 



To generate the CGC initial conditions we use the model from 
gluons are sampled using the kx factorization ansatz 69] 
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681 ] . In that model 



dN n 



„max 

Pt pt 
d 2 pr '' A ~ 



d 2 ki 



a s (f) A Xi 



(p T + k T ) 



(Pr-k 



,(30) 



d 2 r T dy N 2 - 1 J p z T J 4 ~' T "y-" 4 jr~y--~, 4 
where N c = 3 stands for the number of colors, X\p = px exp(±y) / \fs for the light cone 
momentum fractions and \fs for the center of mass energy. 
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According to the KLN (Kharzeev-Levin-Nardi) ansatz 
distribution function <p(x, k 2 -) is related to the saturation scale Q sat via 

1 Q 



7l\ . the unintegrated gluon 



(x, k\\ r T ) ~ 



) 2 

' sat 



a s(QLt) max(Q2 at; k 2 ) 



(31) 



The normalization of Equation (31) and therefore also (30) is determined by the experimen- 



tally measured particle multiplicity at mid-rapidity in central collisions at RHIC. For LHC 
it is normalized in such a way that it is in agreement with the yield from PYTHIA. 
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(a) (b) 
Figure 12. (Color online) Rapidity distribution of parton number dN/dy (a) and their transverse 
energy d-Ey/dy (b) for PYTHIA, mini-jet and CGC initial conditions in a central Au+Au collision 
at RHIC. 




Figure 12 compares the initial particle and energy distributions of all partons for the 
three used models as a function of rapidity. Although the particle number at mid-rapidity 
in all three models is nearly the same, the energy distributions differ considerably. The 
small value from PYTHIA is a consequence of omitting beam remnants like diquarks, since 
BAMPS works on a partonic level. A future project will be to include string fragmentation 
of diquarks and partons in order to transport part of their energy into the partonic phase. 



Figure 13 shows the same distributions for LHC The model for the CGC only samples 



gluons with rapidity smaller than 6, which is, of course, not a physical effect. However, 
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since we are interested in the mid-rapidity region that is not a problem. PYTHIA yields 
also for LHC the smallest initial energy distribution. The particle number distributions are 
comparable with various ot 



rer models 



721 ] . However, there are huge uncertainties due to 



shadowing of 30 — 40 % [62|, |63|, which also affects secondary charm production in the QGP 
due to the unclear initial properties of the medium. On the other hand primary charm 
production in initial hard scatterings is influenced by shadowing as well (cf. next section). 



C. Charm production in hard processes 



The most general expression for the double differential cross section for charm pair pro- 
duction at a collision of the hadrons A and B is given by the partonic cross section for the 
orocess aV,^ C c convoluted with the parton distribution functions (PDF) fa in the hadrons 



39fl: 



da. 



AB 

J7J J? — cc 

c d 3 p c d 3 p 5 



^ / dx 1 dx 2 ff(xi, fi F )ff(x 2 , fi F )E c E, 



ddjj^ccjx^A, x 2 Pb, M c , fi R ) 
d 3 p c d 3 p 5 



(32) 



The indices i and j denote the partons from A and B. fx F and ixr are the factorization 
and renormalization scale, respectively, which are mostly chosen to be equal, fx = fx F = jx R , 



where fx = 2M C 
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73 



74j or fx = + M c 2 [41 



m 



75] are common choices. 



1. Leading order charm production within the mini- jet model 



In leading order charm quarks are produced in the two processes 

9 + 9 -> c + c 

q + q^c + c. (33) 

Experimentally, due to the confinement one cannot measure single charm quarks core 
but only mesons with charm content H(cq), H(cq) or K(cc). The invariant cross section for 
an open charm production process A + B — >■ H + H with two charmed mesons in the final 



state reads 



39] 



da A ^ 

HH 



H 



d?p H d 3 pjj 



fdx,dx 2 dz c dz-M^)^- 



D H /c(Zc) D H/c( Z c 



(34) 



-5 4 (Pl +P2-Pc-P, 
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with 

C( Xl ,x 2 ) = ff{ Xl ) ff{x 2 ) ^£ + [tf ff{x 2 ) + #(sO f q B (x 2 )} . (35) 

Here a/J is the center of mass energy of the partons, which is related to the center of mass 
energy of the hadrons \fs via s = xix 2 s. The partonic cross sections are given in Section HT1 
The fragmentation function Dh/ c { z ) describes the hadronization of the charm quarks, 
with z = \ph\/\Pc\- Since fragmentation only affects the momentum distribution but not the 



total cc cross section 



39j , in which we are interested, the momentum of H is assumed to be 



equal to the charm quark momentum, that is, z — 1 and Dh/ c {z) = 6(1 — z). Thus, taking 



energy conservation into account, Equation (34) can be written as 



dai B 

xix 2 C(x 1 ,x 2 ) , (36) 



dp^dy c dy e 

y being the rapidity. In the center-of-mass frame, the charm and anti-charm quarks hold 
the same transverse momentum p^. 

Integration of the previous equation leads to the total cross section 

If f f da A - B 

o a -- 2 J <1A J iVc J * j^pj^ . (37) 

? ,min ,,min 

The factor 1/2 results from the method of charm number counting, which is widely adopted 
in the literature. The charm quark number corresponding to the cross section should not be 
the total number of particles (charm plus anti-charm quarks) but the number of charm pairs. 



The number of produced charm pairs and the cross section is simply related via 



21, 



44 1 



N c5 = a^T AB (b) . (38) 



The limits of the integrals over the rapidities in Equation (37) can be obtained from the 



definition of the rapidity and some kinematic considerations, 



with xt = 2M T /y/s. 



y™ ax/mia = ±ln(—-e ± y) , (39) 
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From equations (34) and (35) it is obvious that the charm production cross section is 
fundamentally dependent on the PDFs. Therefore, using different PDFs provided by the Les 
Houches Accord Parton Density Function (LHAPDF) group [t6 | we explore this dependency 
in detail, as listed in ITable II 

For that study we also vary the charm mass and the factorization and renormalization 
scale. The calculations are done with a running coupling a s . In contrast, a fixed coupling of 



0.3 yields much smal 



and experimental data 



77 



er results. Additionally, results from PYTHIA (see next section) 
- 79( 1 are shown in the table. The reason for the strong deviation 
between both experiments is currently under intensive investigation. 
Our findings are in good agreement with the literature. In 



the charm production 



cross section in LO is calculated to a 



NN 



shadowing into account. 



160 fib using the HIJING model and taking 
441 ] finds a value between 133 and 153 fib depending on the choice 
of the parton distribution functions. In [80] , the cross section lies between 300 and 750 fib, 
according to whether shadowing and/or strong color electric fields are incorporated. 

Our results in ITable II as well as those from the previously mentioned LO models are 
strongly dependent on various parameters. Mainly, the dependence on the scales shows 
that LO calculations are not sufficient to describe these processes. Indeed, compared to 
experimental data, the LO results are much too small. In contrast, next-to-leading order 
(NLO) calculations are closer to the data [4l|. However, errors due to uncertainties in PDFs, 
scales and charm mass are still significantly high. 

ITable III lists our predictions for the charm production cross section of nucleon-nucleon 
collisions at LHC energy. Again, there is a considerable dependency on the aforementioned 
parameters. 36] calculates the cross section for LHC to a^ 5 N = 5750 fib and [s^] to 6400 fib. 



In 



44j it varies between 2000 and 7000 fib, depending on the PDFs chosen. 



2. Charm production with PYTHIA 

We also used PYTHIA to estimate the number of produced charm quarks in a nucleon- 
nucleon collision. The corresponding cross sections are listed in Tables HI and HTl Applying the 
scaling prescription introduced in Section HV A| the number of initial charm quarks in a heavy 
ion collision can be computed from this value. Since they are produced in hard processes, 
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65 
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66 
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43 
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Table I. LO pQCD cross sections for charm pair production in nucleon-nucleon collisions at y/s 



200 GeV with a running a s , N* = 3, Aqcd = 346 MeV 12j] and a K factor of K = 1 for various 
parton distribution functions, renormalization (/xr) and factorization (hf) scales and charm masses 



M c . Results from PYTHIA and experimental data 
1.5 GeV is used. 



77 



79(] are also listed. In PYTHIA M c 
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Table II. As in ITable II but at LHC energy of yfs = 5.5 TeV. The results for GRV981o with a 
factorization scale of fip = \jp\ + -^c anc ^ f° r PYTHIA cannot be calculated, because these PDFs 
are not designed for such large scales. 

their number scales with the number of binary collisions^] ITable I Hi lists the number of 
initially produced charm quarks in Au+Au collisions at RHIC according to PYTHIA for 
various parton distribution functions. Again, the uncertainties of using different PDFs are 
reflected. The mean value lies at about 9 pairs. 

1 Actually, the scaling factor is a bit smaller if one takes shadowing into account as is done in this calculation 
(cf. Section WM- 
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Number of charm pairs 


CTEQ51 (LO) 


[81] 


8.9 


CTEQ61 (LO) 


[82] 


9.2 


CTEQ6m (MS) 


[82] 


13.6 


MRST2001LO 


[83] 
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MRST2007LOmod 
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HERAPDF01 


[85] 


12.3 


GJR08 (FF LO) 


[86, 87J 


3.0 


GRV98 (LO) 
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Table III. Number of charm pairs produced in primary hard scatterings in central Au+Au colli- 
sions at RHIC for some parton distribution functions by sampling nucleon-nucleon collisions with 
PYTHIA and scaling to Au+Au collisions (cf. Section [IV A j) . 



According to 89| 2 — 6.5 charm quarks are produced in initial hard collisions at RHIC 
in LO, dependent on the PDFs, charm mass and shadowing. {38] estimates this number 
to 2 using the mini-jet model and a phenomenological factor of K — 2. 90[ extrapolates 
from p+p collisions in NLO calculations 8.7 produced pairs. The authors also calculate the 
number of charm quarks at mid-rapidity to about 3 pairs, which is in good agreement with 
our results of about 2 for CTEQ61 and 3 for CTEQ6m. In another NLO calculation 91 1 
the total number of charm quarks varies between 8 and 13. Thus, recent results from the 
literature are in good agreement with our results from PYTHIA. 

Experimental data for the differential charm production cross section at mid-rapidity 
in a nucleon-nucleon collision at yfs = 200 GeV da^ 5 N /dy is available both from STAR 
and PHENIX. The former one measured a value of da^ 5 N / dy = 300 ± 130 fib [79) and 



the latter of 123 ± 47 fib [78!], again showing big deviations. The rapidity distribution 



of the charm production cross section simulated with PYTHIA is depicted in Figure 14 
together with the experimental data points and the pQCD calculations from Section TlV C 11 
The charm distribution from the LO pQCD calculation lies far below the experimental 
data. Hence, higher order corrections should be taken into account or the introduction of a 
phenomenological factor of K > 2 is necessary. The results from PYTHIA agree well with 
the PHENIX data point, which is a bit peculiar, since PYTHIA is based on LO pQCD cross 
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Figure 14. (Color online) Charm production cross section da^ N /dy as a function of rapidity y in 



A and pQCD, respectively, for 



78 



79J]. The pQCD calculation 



a nucleon-nucleon collision at RHIC energy simulated with PYTH 
the PDFs CTEQ61 and CTEQ6m together with experimental data 

is done in LO with fj, F = fj, R = yJp 2 T + M 2 , N f = 3, M c = 1.5 GeV, A QC d = 346 MeV [12] and 
K = 1. 



PDF 


Reference 


Heavy 
Charm 


^uark pairs 
Bottom 


CTEQ61 (LO) 
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62 


7.2 


CTEQ6m (MS) 


[82] 


66 


6.9 


MRST2007LOmod 
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8.9 



Table IV. As in lTable IIIl but for central Pb+Pb collisions at LHC. 



sections. However, PYTHIA is tuned with a running coupling and K factors in order to 
describe experimental data well. 

Although CTEQ6m reproduces the data better than CTEQ61, we will use in the following 



the latter PDF set, because it is designed for LO event generators such as PYTHIA 82]. 
Therefore, we can use PYTHIA also for a reliable sampling of the initial gluon distribution. 
ITable IVl lists the number of charm and bottom quarks produced during initial nucleon 
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Figure 15. (Color online) As in Figure 14 but for LHC. 



scatterings in Pb+Pb collisions at LHC according to PYTHIA. The results from different 
PDFs deviate not as much as for collisions at RHIC. Nevertheless, there are still big system- 
atic errors due to uncertainties in mass, shadowing, factorization and renormalization scale, 
although these are not reflected in the table. 

predicts 34 produced charm pairs in Pb+Pb collisions at LHC within the mini-jet 



model. The NLO value from 



901 ] was estimated to 450, but reduced to 67 — 150 after taking 



newer PDFs and shadowing into account [9JJ. [9l| also makes a prediction for bottom 
quarks, of which about 5 pairs should be produced at LHC. 

The charm production cross sections da^ 5 N /dy as a function of rapidity y for nucleon- 
nucleon collision at LHC energy simulated with PYTHIA and LO pQCD, respectively, with 



different PDFs are shown in Figure 15 It is surprising that the curves from PYTHIA and 
LO pQCD do not differ very much for CTEQ6L. For CTEQ6M, however, the curve from 
PYTHIA is nearly identical to that from CTEQ6L, but it differs by more than a factor of 
3 from the LO pQCD result for the same PDF. From this surprising result, which can also 
be observed in ITable 1V\ one sees again the big uncertainties due to different parameters, 
models, PDFs and shadowing. The NLO calculation including shadowing from [44( lies 
about a factor of 2 above the PYTHIA yield with CTEQ6L. 
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V. HEAVY QUARK PRODUCTION IN THE QGP 



A. Charm production at RHIC 

In this section, charm production during the quark gluon plasma phase in central Au+Au 
collisions at RHIC will be investigated within the framework of the parton cascade BAMPS 
(cf. Section[II])- In contrast to Section IIH] where we limited ourselves to a static medium with 
just two processes, we consider now the full BAMPS simulation of the expanding fireball 
with all interactions from (2) In that framework, because of the gg -H- ggg processes, 



rapid thermalization [34||, for instance, or the build-up of the elliptic flow [92] of the gluonic 
medium can be explained. 

The initial distribution of the charm quarks is simulated using PYTHIA with the parton 
distribution function CTEQ61 and scaled to heavy ion collisions as described in Section 
II V Al Since light quarks are not implemented in BAMPS yet, they are treated as massless 
gluons to take conservation of energy and particle number into account. In this gluonic 
medium charm quarks are produced by gluon fusion g + g — > c + c and interact with gluons 
in elastic scatterings within our model. For both processes, LO cross sections (cf. Section 
HI]) and a constant coupling of a s = 0.3 are used. First, we study charm production without 
multiplying the cross section by any K factor, although we will occasionally employ K = 2 
later. Charm annihilation c + c— > g + g is very unlikely due to the small number of charm 
quarks which are produced and can be neglected [27 ]. 

The primary sources of produced charm quarks at RHIC are initial hard parton scatterings 
during nucleon-nucleon collisions. With PYTHIA and CTEQ61 the initial yield amounts 
to 9.2 charm quark pairs as was shown in the previous section. The time evolution of the 
number of charm quarks during the quark gluon plasma phase within our BAMPS simulation 



is depicted in Figure 16 The number of charm pairs increases slightly in the first 1 — 2 fm/c 
by 0.3 pairs on average due to the high initial temperature of the plasma. Subsequently, 
the number saturates at a value of 9.5 charm pairs after 5 fm/c, which indicates a marginal 
charm production of only 3 % of the total final charm quarks in the QGP. To illustrate this, 
the full range of the y-axis is shown in the plot in the inset. Additionally, the experimental 
value of charm pairs, which is obtained from analyzing D meson decays, is also shown. 
However, we calculated N c£ and its errors from N c£ /Taa and Taa since these are the only 
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Figure 16. (Color online) Number of charm quark pairs produced in a central Au+Au collision at 
RHIC according to BAMPS. The initial parton distribution is obtained with PYTHIA. Note that 
only a small interval of the y-axis is shown. Therefore, we also depict the full range in the inset. 



In addition, the experimental value for the number of final charm pairs is plotted 



93] 



values given in [93]. The number of charm pairs estimated with BAMPS and PYTHIA is 
below the experimental value, but still within the errors. This deviation results from the 
initial charm number obtained with PYTHIA and CTEQ61, which is also smaller than the 



experimental value in proton-proton collisions (cf. Figure 14 ). 

In addition, experimental data show that the number of produced charm quarks at RHIC 
scales with the number of binary collisions 79J, |93j, |94|. That is, charm quarks are mainly 



created in initial hard parton scatterings and not in the QGP, which is in good agreement 
with our findings. 

The small charm yield in the QGP at RHIC confirms our results from Section IIII B 11 
where we estimated a huge chemical equilibration time scale for charm quarks in a gluonic 
box at RHIC temperatures. 

The number of charm quarks produced in the QGP is very sensitive to the initial gluon 
distribution. Therefore, we employed in addition to PYTHIA the color glass condensate 
(CGC) and the mini-jet model for the initial gluon distribution (cf. Section [IV Bp . Both 
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Figure 17. (Color online) As in Figure 16, but with PYTHIA, CGC, and mini-jet initial conditions 
for gluons. For better comparison, the initial charm distribution from PYTHIA is used for all 
models. 



models lead to an increase of a factor of 2.5 in the number of charm pairs produced in the 



QGP compared to PYTHIA, as is shown in |Figure 17] This larger yield is a result of the 
larger initial gluon energy density of both models compared to PYTHIA, as discussed in 
Section HVBl 

Often, the charm production cross section is multiplied by a K factor of 2, which leads 
(naturally) to twice as many charm quarks. Nearly the same effect is attained by altering 



the charm mass from 1.5 GeV to 1.3 GeV, as depicted in Figure 18 



Either changing the initial gluon distribution, choosing a smaller charm mass, or employ- 
ing a K factor can enhance the charm production by about a factor of 2, but the ratio of 
charm quarks produced in the QGP and total charm yield is still just 6 %. A combination of 
all these parameter changes (minijet initial conditions with high initial energy density, small 
charm mass of M = 1.3 GeV and K factor of 2) yields 3.4 charm quarks which are produced 
during the QGP phase. This corresponds to 27 % of the final charm quarks, a value being 
the upper limit of charm production in our model. 

In ITable Vl our results are summarized. In addition to the total charm numbers, we list 
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Figure 18. (Color online) As in Figure 16, but varying the charm mass M and the K factor for 
the cross section of gg — > cc. 



IC 


M c 


K 


dNr*. 
dy 


N c - C 


PYTHIA 


1.5 


1 


1.68 


9.5 


Mini-jets 


1.5 


1 


1.74 


10.0 


CGC 


1.5 


1 


1.74 


9.9 


PYTHIA 


1.5 


2 


1.74 


9.8 


PYTHIA 


1.3 


1 


1.71 


9.8 


Mini-jets 


1.3 


2 


2.16 


12.6 



Table V. Comparison of the total charm pair yield N c e and the number of charm pairs at mid- 
rapidity (y G [—0.5, 0.5]) ^j^ 2 - after the QGP phase at central Au+Au collisions at RHIC simulated 
with BAMPS. Masses M c are in units of GeV. K stands for K factor, and IC for initial conditions. 



the numbers at mid-rapidity. ITable VIl gives an overview about charm production at RHIC 
in other models. The big deviations of the various models result from the influence of many 
parameters, which are often chosen differently in different models. For instance, there is 
disagreement regarding the temperature of the QGP, thermalization time scale, constant or 
running a s , thermal masses of quarks and gluons, volume, energy density, K factor or charm 
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Reference 


^cc.QGP 


M c 


K 


IC 


Comments 


Present work 


0.32 


1.5 


1 


PYTHIA 


Gluonic medium, a s = 0.3 


Present work 


0.80 


1.5 


1 


Mini-jets 


Gluonic medium, a s = 0.3 


Present work 


0.77 


1.5 


1 


CGC 


Gluonic medium, a s = 0.3 


Present work 


0.63 


1.5 


2 


PYTHIA 


Gluonic medium, a s = 0.3 


Present work 


0.61 


1.3 


1 


PYTHIA 


Gluonic medium, a s = 0.3 


Present work 


3.4 


1.3 


2 


Mini-jets 


Gluonic medium, a s = 0.3 


[40] 


0.97 


1.5 


1 


HIJING 


Gluonic medium 


[40] 


1.38 


1.5 


1 


HIJING 


Quarks and gluons with thermal masses 


[40] 


3.2 


1.2 


1 


HIJING 


Gluonic medium 


[40] 


4.9 


1.2 


1 


HIJING 


Quarks and gluons with thermal masses 


[40] 


0.01 


1.5 


1 


Minijets 


Gluonic medium 


[90] 


1 


1.2 


1 


Hydro 




[89] 


3.8 


1.5 


1 


Therm. 


running coupling a s 


[89] 


39 


1.2 


2 


Therm. 


running a s 


[89] 


11 


1.5 


1 


Therm. 


constant a s (M 2 ) = 0.37 


[89] 


120 


1.2 


2 


Therm. 


constant a s (M c 2 ) = 0.42 


[38] 


~ 2 


1.5 


2 


HIJING 




[36] 


~ 0.2 


1.5 


2 


HIJING 





Table VI. Comparison of charm pair yield -/V^qgp in the QGP at RHIC from various models. 
Masses M c are given in GeV. K stands for K factor and IC for initial conditions. 



mass. Within our results we also saw the sensitivity of the production rates concerning some 
of these parameters. 

Another interesting topic we want to address is the question of whether charm quarks are 



chemically equilibrated. For that, we investigate the fugacity defined in Equation (24) for 
charm quarks which are located in central tubes with radius r = 2 fm in transverse direction 
and longitudinal boundaries at space-time rapidity rj = ±0.5. The equilibrium value for 



the charm quark number n^ em ' eq is computed with the effective temperature of the gluonic 



medium. Both variables are shown in Figure 19 for all three initial models. The temperature 
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Figure 19. (Color online) Evolution of the effective gluon temperature (a) and the charm quark 
fugacity (b) in the center of the collision (r < 2fm and r\ E [—0.5, 0.5]) at RHIC for different initial 
models. The insets show the evolution of both quantities for PYTHIA initial conditions until 
8fm when the temperature drops below the phase transition temperature of about 150 — 200 MeV 



95J497] 
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is quite large for RHIC energies, because we consider only a gluonic medium. If we took 
also light quarks into account, the temperature of the medium would be smaller. 

The number of charm quarks is below the equilibrium value and increases with time. As 
we saw in Section IHIB 1[ however, their equilibration time scale is by far too large to see a 
significant rise before hadronization. In addition, the production rate decreases with time 



due to the accompanied temperature decline [50]. If the system evolves for a longer time, 



as is shown in the insets in Figure 19 the fugacity increases dramatically because of the 
decreasing temperature. At the temperature of 161 MeV the fugacity is about 4, which is of 



the same order as the result from 



27, 



98| for the same temperature. However, in 27J, |98| the 



fugacity is computed on a hadronic level from the equilibrium number of charmed hadrons, 
whereas in our approach the fugacity is calculated on a quark level from the equilibrium 
number of charm quarks. Taking this and the difference of the considered volume into 
account, it is understandable that the values of the fugacity are not exactly the same. 

At this low effective temperature of 161 MeV, the energy density in the considered volume 
is about 0.3 GeV/fm 3 (in a totally equilibrated gluonic medium it would be 0.43 GeV/fm 3 ). 
As a note, if the energy density in a cell in BAMPS drops below the critical energy density of 
1 GeV/fm 3 the particles stream freely and do not interact anymore. In the present simulation 
this energy density in the considered volume corresponds to T w 220 MeV (cf. T « 200 MeV 
in an equilibrated medium). 



B. Charm production at LHC 

Due to the higher center of mass energy at LHC than at RHIC, more charm quarks are 
produced during the QGP phase and the yield from initial nucleon-nucleon scatterings is 
also much higher at LHC. However, since the charm production in the QGP increases almost 
exponentially with the QGP temperature but the initial charm yield only logarithmically 
with the collision energy [50 1, the importance of the former rises significantly. 

For the initial gluon distribution we use again PYTHIA, the CGC, and the mini-jet 
model. In Section IIVC2I we estimated with PYTHIA and CTEQ61 about 62 charm pairs 
for the initial yield at LHC, which evolve in the QGP as a function of time as shown in 



Figure 20 With initial conditions from PYTHIA about 11 charm pairs are produced, which 



corresponds to about 15 % of the total final charm quarks. In contrast to RHIC, where 
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Figure 20. (Color online) Number of charm pairs in a central Pb+Pb collision at LHC simulated 
with BAMPS with PYTHIA, CGC, and mini-jet initial conditions for gluons. 



IC 


M c 


K 


dN«r. 

Ay 


N c5 


PYTHIA 


1.5 


1 


7.8 


74 


Mini-jets 


1.5 


1 


14.7 


118 


CGC 


1.5 


1 


10.6 


88 


PYTHIA 


1.5 


2 


9.0 


85 


PYTHIA 


1.3 


1 


8.7 


83 



Table VII. As in lTable V[ but for central Pb+Pb collisions at LHC. 

nearly all charm quarks are produced in initial hard parton scatterings, at LHC a sizable 
fraction is created in the QGP. For initial gluon distributions from the CGC or the mini-jet 
model this fraction is even higher. In the latter framework the charm yield in the QGP is 
actually comparable to the initial yield (47% of the total charm comes from the QGP). The 
introduction of a K factor of 2 or the lowering of the charm mass to M c = 1.3 GeV increases 
these values again by a factor of 2. 

ITable VLTl summarizes our results and ITable VLU1 compares them to other models. As 
for RHIC (cf. previous section), most of these predictions from other models are only 
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Reference 


N C c, QGP 


M c 


if 


IC 


Comments 


Present work 


11 


1.5 


1 


PYTHIA 


Gluonic medium, a s = 0.3 


Present work 


55 


1.5 


1 


Mini-jets 


Gluonic medium, a s = 0.3 


Present work 


26 


1.5 


1 


CGC 


Gluonic medium, a s = 0.3 


Present work 


23 


1.5 


2 


PYTHIA 


Gluonic medium, a s = 0.3 


Present work 


20 


1.3 


1 


PYTHIA 


Gluonic medium, a s = 0.3 


Present work 


38 


1.5 


2 


PYTHIA 


Gluonic medium, a s = 0.3 


[40] 


43 


1.5 


1 


HIJING 


Gluonic medium 


[40] 


94 


1.5 


1 


HIJING 


Quarks and gluons with thermal masses 


[40] 


102 


1.2 


1 


HIJING 


Gluonic medium 


[40] 


245 


1.2 


1 


HIJING 


Quarks and gluons with thermal masses 


[40] 


21 


1.5 


1 


Minijets 


Gluonic medium 


[90] 


23 


1.2 


1 


Hydro. 




[QQl 


O 


1.0 


i 
i 




t = 0.25 fm/c, T = 1.02 GeV, X g = 0.43, 
X q = 0.082 


[50] 


(7) 


1.3 


1 


Therm. 


T = 700 MeV, NLO, at mid-rapidity 


[50] 


(33) 


1.3 


1 


Therm. 


T = 750 MeV, NLO, at mid-rapidity 



Table VIII. As in lTable Vll but for central Pb+Pb collisions at LHC. 



for thermally produced charm quarks, whereas in our calculations the prethermal charm 
production within the first ~ 0.5 fm/c plays a crucial role. 



Figure 21 shows the evolution of the plasma temperature and the charm quark fugacity at 
LHC. The charm quark fugacity being below 1 and the smaller chemical equilibration time 
at this temperature (cf . Section IHIB 2j) are the reasons for the significant charm production 
during the QGP phase at LHC. The inlays show that the fugacity increases dramatically, 
when the temperature approaches its value at the phase transition. At a temperature of 
161 MeV the charm quark fugacity is about 4, which is considerably smaller than the fugacity 
from [27J based on charmed hadrons. 
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Figure 21. (Color online) As in |Fignrc 19 but for LHC. 



C. Bottom production at LHC 



In Section IIVC2I we estimated the initial bottom yield at LHC with PYTHIA and 



CTEQ61 to about 7.2 pairs. Figure 22 shows the evolution of the number of bottom pairs, 
which is only slightly dependent on time. According to our simulation only 0.01 bottom 
pairs are produced on average in the QGP at LHC, which corresponds to about 0.2% of 
the total final bottom quarks. Consequently, bottom production in the QGP phase is neg- 
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Figure 22. (Color online) Number of bottom pairs produced in central Pb+Pb collisions at LHC 
simulated with BAMPS. The initial parton distribution is sampled with PYTHIA. 

ligible. This is a result of t he hu ge bottom mass Mb 3> T, which is much larger than the 



temperature of the medium 



100]. 



Therefore, bottom quarks are a promising probe at the LHC. On the one hand their 
number is large enough to be measurable; on the other hand one knows that bottom quarks 
are produced initially. Due to this information one will be able to draw conclusions about 
the early phase of the QGP from observables like the elliptic flow or the nuclear modification 
factor. 



VI. CONCLUSIONS 

We have studied charm and bottom production in heavy ion collisions at RHIC and LHC 
using PYTHIA and the mini-jet model for the primary yield and the partonic transport 
model BAMPS for secondary production in the QGP. For that, we implemented the LO 
2 — > 2 processes for heavy quark production, their annihilation and their interaction with 
gluons. The charm production in the QGP at RHIC can be neglected, but at LHC we 
expect that open charm will no longer scale with the number of binary collisions due to a 
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large charm production in the QGP. 

The BAMPS simulation of charm production in a static gluonic medium was compared 
with the analytic solution of the corresponding rate equation, where an excellent agreement 
was found. From these calculations we estimated the order of magnitude of the chemical 
equilibration time scale for charm production. For a medium with an initial temperature 
of 400 MeV - the approximate temperature of the QGP at RHIC - we found a time scale 
of about 700fm/c, which is much larger than the lifetime of the QGP hinting of a small 
charm production at RHIC. For an LHC temperature of 800 MeV, we estimated the chem- 
ical equilibration time scale to be 30fm/c. For slightly higher temperatures, it is even 
comparable to the lifetime of the QGP at LHC. Consequently, one can expect a significant 
charm production in the QGP at LHC, which is, however, still strongly dependent on the 
initial conditions. 

The production of heavy quarks in initial hard parton scatterings during nucleon-nucleon 
collisions was estimated with PYTHIA and the mini-jet model in LO pQCD, revealing a 
strong sensitivity to the chosen parton distribution functions, heavy quark mass, factoriza- 
tion, and renormalization scale. Charm production within the mini-jet model lay about a 
factor of 5 below the PHENIX data, while PYTHIA's value was close to the experimental 
data. 

Based on this initial heavy quark yield, we studied the production of charm quarks in 
central heavy ion collisions at RHIC and LHC and bottom quarks at LHC in a full space-time 
simulation of the QGP with BAMPS. In addition to the 9.2 initially produced charm pairs 
with PYTHIA and CTEQ61, between 0.3 and 3.4 pairs were created in the QGP, depending 
on the initial gluon distribution (from PYTHIA, CGC, or mini-jets), charm mass, and K 
factor. Because of this small production, the charm quark number during the QGP phase 
is below its equilibrium value and only rises above this value when the temperature of the 
medium drops below 200 MeV. 

For Pb+Pb collisions at LHC, the number of initially produced charm pairs is 62 according 
to PYTHIA with CTEQ61. Subsequently, in the QGP, the yield is significant: between 11 
and 55 pairs are created, depending on the model for the initial gluon distribution. For the 
LHC, we estimated the initial bottom yield as 7.2 pairs, whereas their production in the 
QGP can be neglected. 

In conclusion, heavy quarks are an interesting probe for the early stage of heavy ion 
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collisions due to their strong sensitivity on the initial conditions and the properties of the 
QGP within the first fm/c. However, it is crucial to disentangle the origins of heavy quarks: 
whether they were produced in primary or in secondary production. If one knows the ratio 
of both, the measured total yield can reveal information about important variables in the 
QGP: The primary yield is strongly influenced by the parton distribution functions, mass, 
scales and also shadowing; the secondary production depends on the temperature of the 
medium, which is again dependent on the initial conditions including shadowing. Therefore, 
both yields are correlated and can be disentangled by comparing the effects of using different 
models for initial conditions and shadowing on the total heavy quark yield within an accurate 
description of the QGP. In the present article, we have studied the impact of different initial 
conditions but considered shadowing just by lowering the number of binary collisions. In 
a future investigation, we will take more sophisticated models for shadowing into account, 
which will also affect primary heavy quark production. 



Other observables, that can give great insight into the dynamical properties o 
are the e lliptic flow and the nuclear modification factor Raa of heavy quarks 



theQGP, 



101 



20-23. 



1021 ] . Experimentally, a strong elliptic flow of heavy flavor electrons and a small Raa 



have been observed 



19], which could not be described wit 



pQCD 103| . However, next-to-leading order calculations 



lin t he leading order regime of 



23 



10 lfl that also include radiative 



energy loss are very promising for explaining the data. A future task will be to investigate 



the energy loss and these observables with our model, implementing also li ght quarks, 



order corrections, a running coupling, and an improved Debye screening 



24 



25 



ligher 



104 



105]. 



For that, BAMPS is a well suited model, since the framework for 2 <H- 3 interactions is 
already imple mented |33l| . Extending recent BAMPS studies concerning v 2 and Raa of 
gluons [ill . I92I 106 , 107| | to the heavy flavor sector will reveal further information about 



interactions in the QGP. Moreover, we want to study the effect of hadronization on the 
elliptic flow or Raa, which could b e done, fo r ins tance, in a fragmentation, coalescence or 



statistical hadronization model 



27 



103 



108 



. In this framework J/ip creation can be 



investigated as well. If n ecessary, we will also consider other schemes such as AdS/CFT 



correspondence 



110 



llll ] in addition to pQCD. 
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